Probing the interior of neutron stars with gravitational waves 
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We show here how the internal structure of a neutron star can be inferred from its gravitational 
wave spectrum. Under the premise that the frequencies and damping rates of a few w-mode oscilla- 
tions are found, we apply an inversion scheme to determine its mass, radius and density distribution. 
In addition, accurate equation of state of nuclear matter can also be determined. 
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INTRODUCTION 

Despite that direct detection of gravitational waves 
from violent astronomical events such as coalescing bi- 
nary neutron stars (or black holes) and supernovas has 
not yet been realized, it is generally believed that it can 
be done in a few years. Gravitational wave detectors of 
different designs, including resonant antennas (e.g. EX- 
PLOPER and NIBOE), LIGO and LISA (currently under 
construction) , complement the sensitivity of one another 
and joint observations among them are likely to increase 
confidence and accuracy in detection (see, e.g. 0,0] and 
references therein). In addition to providing cogent sup- 
port for general relativity, gravitational wave detection 
indeed opens up a new channel for us to survey stellar 
objects far away from our galaxy Q,0|. 

Neutron stars and black holes are undoubtedly two 
major sources of gravitational waves. However, unlike 
black holes that are completely specified by a few physical 
quantities, the structure of neutron stars is still masked 
because the behavior of nuclear matter at high densities 
is not yet completely known. In fact, there are currently 
numerous theoretical equations of state (EOSs) of nu- 
clear matter, leading to different models of neutron stars 
HUH El U, H El To illustrate the differences in these 
EOSs, we show a few of them in Fig.^ including models 
A and C models AU and UT 8], models APR! 

and APR2 

As gravitational waves might carry the imprint of the 
internal structure of a neutron star from which they are 
emitted, a close examination of relevant signals, com- 
monly analyzed in terms of quasi-normal modes (QNMs) 
characterized by complex frequencies a; = u> T — koi |ldllll| . 
could possibly unveil its structure. In fact, the term 
"gravitational wave asteroseismology" has been coined 
to indicate studies along this direction HI [T^. Noting 
the universal behavior in the QNMs of the fundamental 
fluid /-mode and the first polar w-mode, Andersson and 
Kokkotas LL2| showed that the radius and the mass of a 



neutron star can be roughly estimated from simultane- 
ous detection of such wave modes. In addition, it was 
argued that the frequency of the leading p-mode (or the 
axial w-mode) oscillation could be used to identify an 
appropriate EOS describing the neutron star |12.ll3j] . 

In a series of recent papers fl4lll5j ]. we have pointed out 
the physical origin of the universality in the frequencies 
of the /-mode and the w-mode, and worked out a scaled- 
coordinate logarithmic perturbative theory (SCLPT) to 
locate axial w-mode oscillations. The main objective of 
this Letter is to work out an inversion scheme to infer the 
internal structure of a neutron star from its gravitational 
wave spectrum. Under the premise that the frequency u r 
and the damping rate ui; of a few w-mode oscillations of a 
neutron star can be identified from gravitational wave ob- 
servation, we iteratively invert SCLPT to determine the 
mass M, the radius R and the density distribution p(r) 
of the star. In particular, accurate EOS of nuclear mat- 
ter can also be obtained. While the idea underlying such 
inversion scheme reported here is completely generic, we 
will consider and use the spectrum of axial w-mode oscil- 
lations to illustrate our method. Unless otherwise stated, 
geometric units are used in this Letter. 

We first review the universality in QNMs of neutron 
stars 14] . As shown in Fig. [3 where Moj r and Mlo\ 
are plotted against the compactness C = M/R for non- 
rotating neutron stars described by different EOSs, the 
frequency of the w-mode (or /-mode) oscillation approx- 
imately satisfies the following scaling law [13, Il3l Il4 , Il5| : 



Muo = aC 2 + bC + c. 



(f) 



with a, 6, and c being complex constants determined 
from curve fitting. Such universality originates from the 
fact that the mass distribution inside physical neutron 
stars can be approximated by the Tolman VII model 
(TVIIM) 0, whose mass distribution function m(r) = 
4-7T J Q x 2 p(x)dx is given by: 



m c (r) = M 



(2) 
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QNMs of TVIIM manifestly depend only on the compact- 
ness of the star and reproduce the scaling behavior men- 
tioned above 0] . As shown in Fig. TVIIM provides a 
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good approximation to stars with varying EOS and the 
universality in Q can be captured by the best quadratic 
fit to the QNMs of TVIIM, with a = -4.4 - 6.3i, 
b = 3.1 + 1.9i, and c = -0.072 + 0.098i 0. 

As the first step of the inversion scheme, we study the 
frequency of the leading (i.e. the least-damped) axial 



-mode of TVIIM, 



;jy — ic*>xi , and find that 
the ratio w', is in fact a monotonically increasing 

function of M/R Q. Hence, once t^ c) for TVIIM is 

(c) (c) 

known, M/R can be obtained from the ratio dj\ r /u)^ 
and in turn M and R can be found from or Fig. [21 

As TVIIM indeed provides a benchmark for other re- 
alistic stars, we expect that the frequency of the least- 
damped axial w-mode emitted from a realistic neuron 
star, u>i — <jji r — iwn, is close to that of the TVIIM star 
with the same mass and the same radius. Therefore, we 
could go through the procedure outlined outlined above 



with 



.(<=) 



replaced by lj\ to obtain estimates of M, R, 



m(r) and hence p(r) for the star in consideration. 

To gauge the accuracy of this scheme, we show in 
Fig. EJa) p(r) of an APR1 star with C = 0.28. The re- 
sult obtained from the above step (empty circles) agrees 
nicely with the exact numerical data (the solid line). The 
pressure distribution P(r) then follows directly from the 
Tolman-Oppenheimer-Volkoff (TOV) equations 0, UtI ] 
and in turn the EOS P(p) can be found. In Fig. Efb) 
we show the result of such scheme (empty circles) , which 
satisfactorily reproduces the EOS (solid line) except at 
high densities. In addition, as TVIIM is exactly solvable 
the EOS can be found analytically: 
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with po = p(r = 0) = 15C/(87ri? 2 ), 

/(C) = 3C - 3C 2 (1 - p/p ) (2 + 3p/p ) , 
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It is worthwhile to note that a proposal has been put 
forward to invert the EOS from given values of M and R 
[Tsj l . In our scheme we make use of the values of u>i T and 
wii , which are considered as data obtained from gravita- 
tional wave astronomy, to determine M and R, and in 
turn obtain the EOS. Moreover, our scheme can yield an 
analytic expression for the EOS while the proposal in 0] 
has to resort to iterative numerical solution. 

In the following we will further consider the frequen- 
cies of non-leading QNMs, u> q = uj ql — iuj q i (q = 2, 3, . . .), 
ordered in increasing u) qv , and show that they can lead 



to much improved estimates of the mass, the radius, the 
density and EOS. As an incentive to the readers, we first 
display the accomplishment of our scheme. In Table I 
we tabulate the exact values of M and R and those ob- 
tained from our scheme using one or two QNMs for stars 
constructed with different EOSs. The result is truly en- 
couraging, especially for the mass. Hence, we conclude 
that the frequencies of two leading w-mode oscillations 
can readily lead to accurate determination of M and R. 
In Fig. 1 we compare various EOSs and demonstrate that 
it suffices to use only two frequencies (lui and L02) to ac- 
curately reproduce and distinguish these EOSs. Figure 3 
shows p(r) and P{p) obtained from the inversion scheme 
using one, two and three QNMs, respectively. It is clearly 
seen that the result obtained from a combination of u>i 
and u>2 (grey circles) are readily improved, while using 
all three frequencies (dark circles) indeed yields a perfect 
match with the exact values. We have also used uj\ and 
the frequency of a u>n-mode, Wn, to infer the physical 
quantities. However, the result is slightly less accurate. 

The technical details of our scheme are as follows. To 
improve on TVIIM, we write m(r) as: 



m(r) « mo(r; M, R) + pjrrij(r; M, R), (6) 



where 



mj(r;M, R) 



3=1 



15Mr / r \ 2j 



(7) 
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p > 2, and pj (j = 1, 2, . . . ,p — 2) are adjustable param- 
eters to provide the best fit to m(r). As M and R are 
not yet exactly known, they are also considered as free 
parameters. There are consequently p free parameters to 
be determined from the QNM frequencies. As mentioned 
above, m(r) w m c (r) and therefore we will use m c (r) as 
an initial guess for m n (r). 

In the first step of our scheme, approximate values 
of the mass and the radius of a star have been found 
from u>x, which are equal to Mo and Ro, respectively. 
Much improved values of these quantities, M = AIq + Mi 
and R = Rq + R\ , and the mass distribution inside the 
star are now determined from higher-order QNMs as fol- 
lows. Suppose n (n > 1) QNM frequencies of a neu- 
tron star are known. We look for optimal values of pj 
(j = 1, 2, ... ,p— 2), Mi = Pp-\, and R\ = p p that yields 
frequencies with minimal deviations from the measured 
ones. If all of these p free parameters vanish, the cor- 
responding frequencies are uj q °\ Based on SCLPT, the 
first-order change in uiq (q = 1, 2, . . . , n) is 



P a P 

EOLUg \ - 



which is close to the exact change A q = oj, 
plicit expression for c q j can be found in 
expressible as a bilinear form of the wave function. 



15], 



(8) 
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By minimizing the total square-deviation x 2 : 

n 

X 2 ^|A 9 - W «| 2 , (9) 

9=1 

we can obtain a set of p linear equations for fij : 

p 

b k =J2 A kj»j, (10) 

3=1 

where k = 1, 2, . . . ,p , 

n 

b k = ^ (AgC qk + K q C qk ) , (11) 

9=1 

n 

A kj = ( C lkCq + C qk C qj ) . (12) 

9=1 

and we have used an overbar to denote complex conju- 
gation. Consequently, {/ij} can be computed. 

Once {fij} are obtained, we can construct a new model 
star by incorporating these perturbations into mo(r) and 
look for a better approximation of the star by solving (|1(J|) 
again. Through this iterative scheme, we can generate 
rapidly convergent sequences for M, R, p(r) and P(p). 

To demonstrate the robustness of the inversion scheme, 
we have also applied our scheme to infer the EOS of 
a strange star H 0- The EOS of quark matter, P = 
(p — 4B)/3, is given by the MIT bag model 0, where 
B = 57McVfm~ 3 is the bag constant. The crust of 
the strange star, where p < neutron drip point, is con- 
structed with the Baym-Pethick-Sutherland EOS [10. 
The EOS for quark matter is rather exotic and the den- 
sity is even discontinuous across the inner boundary of 
the crust. However, as shown in Fig. the inversion 
scheme still works nicely when four QNMs are used. 

Despite the mass of some neutron star binaries can 
be inferred from their orbital periods, as yet there is no 
generic method to determine their radii. In this Letter 
we have proposed a robust inversion scheme to determine 
the mass, the radius, the mass distribution and the EOS 
of a neutron star from its gravitational wave spectra. We 
expect that our scheme could operate in conjunction with 
gravitational wave detectors in the near future to probe 
the interior of neutron stars. 
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grant (Project ID: 2060260) from The Chinese University 
of Hong Kong. 
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IX I Kill J 


APR1 


4.460/4.475/4.460 


11.83/12.12/11.70 


APR2 


4.173/4.184/4.173 


11.07/11.26/10.91 


AU 


3.843/3.862/3.843 


10.19/10.61/10.08 


A 


3.289/3.282/3.289 


8.723/8.529/8.564 


C 


2.791/2.788/2.791 


7.404/7.185/7.216 


UT 


3.654/3.651/3.655 


9.691/9.745/9.664 



TABLE I: The mass M and the radius R of neutron stars 
constructed with different EOSs are compared with the cor- 
responding values obtained from inversion using one or two 
QNMs. In each entry the values are obtained from TOV 
equation/inversion using u>\ /inversion using u>\ and UJ2 respec- 
tively. 



5 




FIG. 1: Six different EOSs are shown by lines of different 
formats. Grey symbols are the corresponding values obtained 
from inversion scheme using wi and u>2 for stars with C = 0.28. 
Here P and p are in cgs units. 
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FIG. 2: The real and imaginary parts of Mui for the least- 
damped axial ui-mode of six realistic stars (unfilled symbols) 
and TVIIM (stars) are shown as a function of M/R in panels 
(a) and (b) respectively. The dotted line represents the best 
quadratic fit to those of TVIIM. 
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FIG. 3: Panels (a) and (b) depict p(r) and P(p) respectively 
for an APR1 star with C = 0.28. The solid line is the theoret- 
ical value, while the unfilled/grey /dark circles represent the 
results obtained from inversion scheme using one/two/three 
leading axial w-modes. The result obtained from inversion 
scheme using u)i and the frequency of a wn-mode is shown by 
the crosses. 




FIG. 4: Panels (a) and (b) depict p(r) and P(p) respectively 
for a strange star with C = 0.2. The solid line is the theoretical 
value, while the grey triangles/dark triangles/crosses respec- 
tively represent the results obtained from inversion scheme 
using two, three and four leading axial w-modes. 



